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We investigate the dynamics of systems of many coupled phase oscillators with het- 
erogeneous frequencies. We suppose that the oscillators occur in M groups. Each 
oscillator is connected to other oscillators in its group with "attractive" coupling, such 
that the coupling promotes synchronization within the group. The coupling between 
oscillators in different groups is "repulsive"; i.e., their oscillation phases repel. To 
address this problem, we reduce the governing equations to a lower-dimensional form 
via the ansatz of Ott and AntonseiP. We first consider the symmetric case where all 
group parameters are the same, and the attractive and repulsive coupling are also the 
same for each of the M groups. We find a manifold £ of neutrally stable equilibria, 
and we show that all other equilibria are unstable. For M > 3, £ has dimension 
M — 2, and for M — 2 it has dimension 1. To address the general asymmetric case, 
we then introduce small deviations from symmetry in the group and coupling param- 
eters. Doing a slow/fast timescale analysis, we obtain slow time evolution equations 
for the motion of the M groups on the manifold £. We use these equations to study 
the dynamics of the groups and compare the results with numerical simulations. 
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The dynamics of hierarchically structured networks is particularly relevant to 
systems with very large numbers of dynamical units. Examples can be found in 
neuroscience and in the study of social dynamics in large populations. At the 
coarsest level, one can often think of such systems as being based on interactions 
between distinct groups or communities of dynamical network nodes. Thus, it is 
of great interest to investigate the kinds of dynamical behaviors that can result 
when network systems have interacting communities, and to develop techniques 
that may be useful for studying such systems. In this paper we address these 
issues for the illustrative, paradigmatic case in which the nodal dynamics is 
describable within the phase oscillator model. 

I. INTRODUCTION 

It is often observed that real-world networks have groups of nodes that are associated 
strongly with each other and less strongly with nodes in other groups^. Considering situ- 
ations in which the links on such networks represent the flow of signals between interacting 
nodal dynamical units, this type of group structure may fundamentally influence the func- 
tioning, overall performance, and emergent dynamical properties of such systems. In this 
paper we investigate the dynamics of a simple model for the interaction of groups of many 
oscillators. Following Kuramoto, there is an oscillator at each network node whose state is 
given purely by its oscillation phase angle®®. The oscillators occur in M groups and have 
different natural oscillation frequencies. Every oscillator is uniformly and "attractively" 
coupled to all other oscillators in its group, but is also "repulsively" coupled to all oscil- 
lators outside its group. By an "attractive" ("repulsive") network coupling, we mean that 
the coupling promotes (works against) phase synchronization of the two oscillators that are 
linked by the coupling. The general problem of interacting network oscillator groups has 
also been considered^, and the nonlinear dynamics of such groups has been studied for 
the case of groups whose average natural frequencies are not in resonance^ . Our focus on 
within-group attraction and between-group repulsion is partly motivated by our expectation 
that this will provide a particularly clear illustration of the effect of network group structure 
on system dynamics. If we make a loose analogy between the phase oscillator model and 
opinion dynamics, as has been done in the paslP^, we can think of the interactions in the 
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context of a multi-party political system, with within-group attractions representing the 
tendencies for individuals to align with those in the same political party and between-group 
repulsions reflecting the desire of individuals to differentiate themselves from members of 
other parties. 

We consider our model using direct numerical simulations and by the low-dimensional 
reduction (Sec. [n]) of Ott and AntonserP^, which has proven useful for the treatment of a 
wide variety of problems involving the interactions of a large number of phase oscillators^. 



We first consider a symmetric situation (Sec. Ill), in which all parameters describing the 
M oscillator groups are identical, and the attractive and repulsive couplings within and 
between the groups is uniform. It is found, both numerically and analytically, that, above 
a threshold in the coupling, the oscillations within groups may become coherent!^, and the 
groups interact to achieve equilibrium configurations. We also analytically determine all 
possible system equilibria and find that, when M > 3, there is an (M — 2) dimensional (one 
dimensional for M — 2) manifold C of possible neutrally stable equilibria, and that all other 
equilibria are unstable. In order to provide insight into the general asymmetric case (Sec. 



IV), we introduce small asymmetric deviations to the group and coupling parameters. We 



then analyze the situation using multiple-timescale asymptotic^. This results in a set of 
equations describing the slow-timescale evolution of the M group order parameters as they 
temporally evolve and interact on the manifold C of neutrally stable symmetric solutions. 
We also report numerical results testing the applicability of our results. Further discussion 
and conclusion are given in Sec. |Vl 



II. LOW DIMENSIONAL FORMULATION 

In Kuramoto's original formulation, the dynamics of the individual oscillators are given 



by 

^r =t * + £ X>[*;(*)-<W)], (1) 



3=1 

where Q{ is the phase of the zth oscillator [i = 1,2, ...,N), Ui is the zth oscillator's natural 
frequency, and K measures the strength of coupling between oscillators. We consider a 
modified version of this system in which the oscillators are placed into M communities 
(groups) of equal size . Instead of a single coupling constant, we consider a matrix in which 



the aa'th element represents the coupling of each oscillator in group o [a = 1, 2, M) to 
each oscillator in group a'. In this model, the dynamics of the oscillators are governed by 

Aao M t\/tts N/M 

f =«r + £^E°K»r-«r>. (2) 

a' = l j = l 

where d° is the phase of the zth oscillator in group a. The natural frequencies oj° are taken 
from Lorentzian distributions given by 

™ = (M -t)' + A. ' < 3 > 
where A CT and u a are the width and center, respectively, of the Lorentzian distribution for 
group a. 

Just as the dynamics of the original Kuramoto model can be captured by a single complex 
order parameter, the modified system can be understood through the behavior of the M 
group order parameters 

M N/M 

«. = ^E ex P^)- (4) 

i=i 

As shown by Ott and AntonserP^, the group order parameters are attracted to a manifold 
on which they evolve according to the ordinary differential equations (ODEs) 

M 



dotr. 



dt y u u ' u 2 

o'=l 

where the * denotes complex conjugation. Defining 

M 



(A CT - iUJ a )a a - - ^ K <rA a *a> a l - a <r']i ( 5 ) 



Rg = y^-Kp-g/ag/, (6) 



<t'=1 

we rewrite Eq. ^ as 

^ = -(A. - iu a )a a - \{alRl - R a ). (7) 

III. IDENTICAL GROUPS 
A. Equilibria 

We first consider the simplest case, in which = 1 and w ff = Co for all a. In order to 
investigate the effects of repulsion between groups of oscillators, we define the entries of the 



K„ 



(8) 



group coupling matrix to be 

( K :a = a' 
-Kb :a^a' 

where the constant K > represents the overall coupling strength, and b > quantifies the 
degree of repulsion between oscillators in different groups. The average frequency Co can be 
transformed to zero by the change of variables 9{ — > Q{ — uit, and we henceforth set u> = 0. 
[In Section IV we will consider deviations from A a = 1, co a = and Eq. p]).] Defining 



M 



cr=l 



Eq. ^ can be expressed as 
and Eq. ^ becomes 



R a = K(l + b)a a - KbS, 



K 



-otc 



[a 2 a ((1 + b)a* a - bS*) - ((1 + b)a a - bS)] . 



(9) 



(10) 



(11) 



dt "" 2 

In order to characterize possible equilibria of the system, we set da a /dt = 0. Writing a a 
and 5* in polar form as r a e l ^ a and SqC 1 ^ , respectively, we substitute these quantities into 



Eq. (11) and rearrange terms to yield 

K 







r a - - [r CT (l + b)(rl - 1) - rlbS^^ + bS^^} . 



Taking the imaginary part of both sides yields 

A' 



= -65 (^ + l)sin(^-*). 



(12) 



(13) 



This equation gives restrictions on the possible steady-state equilibria in our model. Assum- 



ing K ^ and b ^ 0, Eq. (13) is satisfied if and only if either So — 0, or ijj a — \& = kn for 



some integer k. The latter condition must be satisfied for a = 1, 2, M, implying that, for 
an equilibrium with S ^ 0, the complex order parameters a a all lie on a single line through 
the origin of the complex a-plane. 

B. Equilibria with S = 



We now turn our attention to the equilibria with S = 0. In this case, Eq. (12) becomes 

r„ 1 - —(I +b)(ri - I) . (14) 



l + -(l + 6)(r»-l) 
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Solving for r a yields either r a = or 

2 



Ta = To 



1 - 



1/2 

(15) 



K(l + b)_ 

In Kuramoto's original model, there existed a critical coupling value K c that marked the 



onset of synchronization in the system. For our case, Eq. (15) yields an analogous value, 



K c= — b . (16) 



For K < K C1 the oscillators behave incoherently [r a = is the only solution of Eq. (14)] and 



exhibit no collective behavior. From now on we assume that K > K c , so that the oscillators 
in each group can synchronize with one another, and the groups act as coherent entities. 

For K > K c , an equilibrium with S = can have some group order parameters with length 
zero and some with length r$. It is convenient to visualize the nonzero order parameters 
as lying on a circle of radius r$ in the complex a-plane. Assume that there are m groups 
with r a = and M — m groups with r a = r . If we conceptualize the non-zero order 
parameters as vectors of length r , we can place them tip-to-tail to form an equilateral 
polygon in the plane. Therefore, given M — m groups with r a = r , the set of possible S = 
equilibria corresponds to the set of equilateral (M — m)-gons. For (M — m) — 2, the two 
r a = r groups must have group phases that are diametrically opposed (i.e., separated by 
7r). If (M — m) — 3, the phases of the order parameters must be equally spaced by 27r/3, 
corresponding to an equilateral triangle in the tip-to-tail representation. If (M — m) = 4, the 
tip-to-tail configuration of order parameter vectors corresponds to the sides of a rhombus. 
Since a rhombus has two pairs of parallel sides, there must be two pairs of diametrically 
opposed order parameters, but the angle between them is arbitrary (corresponding to the 
one parameter family of all possible rhombus shapes). For (M — m) = 5 and larger, there is 
more freedom in the specification of equilateral (M — m)-gons (see Fig. [T]), and the pattern 
is not as easily specified as for (M — m) = 2, 3 and 4. In general, assuming (M — m) > 3, 
the set of all possible (M — m)-gon shapes, (and hence the family of all possible equilibria) 
is (M — m — 2) dimensional for (M — m) > 3. That is, these equilibria lie on an (M — m — 2)- 
dimensional manifold. We denote this manifold by C. The case (M — m) = 2 is special; 
in this case, £ has dimension one. For (M — m) = 2 or 3, the dimension of £ is one, 
corresponding to invariance of these equilibria under a rigid rotation of the phases of a a . 





FIG. 1. An S = equilibrium with five groups, and the corresponding equilateral pentagon. \a a \ — r$ < 
l;cr = 1,2,3,4,5. 

C. Equilibria with S ^ 



According to Eq. (13), any solution that does not have S = must have all of the order 



parameters on a single line through the origin. Without loss of generality, we can rotate 
the configuration so that the order parameters all lie on the real line, and the value of 5* is 
positive. Thus we can let ip a = \& = for all a, if we allow r a to take on both positive and 



negative real values. Then Eq. ( 12 ) becomes 

K 



= (1 + b)rl - jbSrl + (l - |(1 + 6) J r a + *bS. 



(17) 



The solutions to this cubic equation give three possible values for r a . We note that since 
S > 0, the product of the roots of this equation, equal to —bS/(l + 6), is negative. If the 
cubic has two complex roots, the third root must be negative to satisfy this condition. Then, 
since the order parameters lie on the real line by assumption, all the order parameters must 
be at the negative real root. But S cannot be positive if all the order parameters have a 
negative value, so we reach a contradiction. Therefore, we conclude that all of the roots of 
this equation must be real. Furthermore, exactly one of the roots must be negative. 



D. Numerical Simulations 



In order to test the validity of these results, we integrated Eq. (11) numerically for each 
oscillator group. We chose values of K and b such that K > K c , and examined the behavior 
of the order parameters for different numbers of groups. We observed that, given arbitrary 




initial conditions for the group order parameters, the system tended to a approach a steady- 
state equilibrium. All typical initial conditions yielded evolutions that tended toward S = 
equilibria in which the lengths of all the order parameters were equal. If we chose the initial 



order parameters so that they lay on a line through the origin and satisfied Eq. (17), we 
obtained stationary states with 5^0. If some of the initial order parameters were zero and 
the initial configuration satisfied S = 0, then we also obtained stationary states. However, 
these configurations were seen to be unstable; any small changes in the order parameters 
caused the system to drift away and evolve towards a state with S = and all r a = r . This 
implies that the only stable states of the system are those in which S — and all r a = r . 



In Sec. HIE, we demonstrate this analytically. 

We also conducted simulations of the full collection of oscillators by integrating Eq. §2§ 
for each oscillator, using a group size of 2000 oscillators. When K > K c , we observed that 
the oscillators resolved into coherent groups with non-zero order parameters of nearly equal 
magnitude. The group phases tended to arrange themselves such that S ~ 0. However, 
because the number of oscillators was finite and their natural frequencies were randomly 
distributed, the average natural frequencies of the groups were not all exactly the same. As 
a result of this asymmetry, the system did not reach a steady-state equilibrium. Instead, 
the group order parameters rotated slowly in the complex plane at varying speeds. The 
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configuration changed continuously, but at every instant of time the system resembled an 
S = equilibrium with all r a = r . As discussed above, the S — condition restricts 
the possible configurations of the system. For M = 2 or 3, the phases of all group order 
parameters changed at the same speed, in order to keep the phases evenly spaced on the 
unit circle. For M — 4, the groups formed two pairs whose phases changed at different 
speeds. For M > 5, the group phases evolved in a more complicated manner, reflecting the 
multitude of possible S = configurations for large numbers of groups. An example with 
M = 3 is given in Fig. [3] 

The absence of stationary states in these simulations is due to the asymmetries that nat- 
urally arise in finite collections of oscillators. Our low-dimensional analysis of the identical- 



groups case does not account for these asymmetries. In Section |IV| we will explain this 
behavior by examining the low-dimensional system in the case in which the groups do not 
have identical properties. 




FIG. 3. A sampling of oscillator phases from a numerical simulation with M = 3. Individual oscillator 
phases 6f appear as solid dots on the unit circle, and the colors (red, green and blue) correspond to different 
oscillator groups. Group order parameters a a appear as larger open circles. 
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E. Stability of Equilibria 



1. Equilibria with S = and r a ^ for all a 

To analyze the stability of the S = equilibria with r a ^ for all a, we introduce a small 
complex perturbation to each order parameter and determine whether the perturbations 



grow or shrink with time. We replace a a with roe^ CT (l + 5a a ) in Eq. (11) and replace S 
with So + SS, where 5S = J2t=i r o^"5a CT . Letting So = and collecting terms that are 
first order in 5a a , we obtain 



d(5a a ) 



—5a„ — 



K 



1 + b)(r 2 5a* a + 2r 2 5a a ) - br e li '°5S* - (1 + b)5a a + —e^'SS 



dt 2 
If we define the vector quantities 

5a = 



iipo 



r 



(18) 



( 5a x \ 

Sao 



y 5a M j 



and £ 



V / 



(19) 



then we can cast the M ODEs in Eq. (18) into the single vector equation 
2 d(5a) 



(r 2 (l + 6))l + 6(cf 



5a — rl 



[1 + 6)1-6 ft 



FT* 



5a 



(20) 



where 1 is the M x M identity matrix, T denotes transposition, and we have made use of 



Eq. (15) to simplify the matrix multiplying 5a. The conjugate of Eq. (20) is valid as well; 



those two equations can be combined into a single vector equation, 

da 
dt 



Aa, 



(21) 



where 




A 




r 2 (l + 6)1 + 6 (f^ 
1 + 6)1-6^ 



[1 + 6)1 -6 (ft 
^(1 + 6)1 + 6^ 



(22) 



We assume an exponential solution a ~ expAt, so that da/dt = Aa for some scalar A. 



Equation (21) is then satisfied if A is an eigenvalue of the matrix A. We can analyze the 



behavior of the system near equilibrium by finding the 2M eigenvalues of A and determining 
whether their real parts are positive, negative or zero. 
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' in \ 


Li 


and a_ = 








I -ml 



We begin searching for these eigenvalues by finding real vectors u satisfying = 0. For 
such a vector w, we see that both 



(23) 



are eigenvectors of matrix A with eigenvalues —[K(l + b) — 2] and 0, respectively. The 
first of these eigenvalues is negative provided that K > K c , which we have assumed. The 
eigenvalue corresponds to neutrally stable perturbations; if a is parallel to an eigenvector 
with eigenvalue 0, then the perturbation will neither grow nor decay in time. 

Noting that 5S can be rewritten as r ^ r 5a, we see that the above condition on u is 
equivalent to requiring SS = 0. Further, since u has real components, the order parameters 
r , oe*^ CT (l + 5a a ) are such that only their magnitudes are perturbed in the case of a+, and 



only their phases are perturbed in the case of a_, where a+ and a_ are as given in (23). 
Because the eigenvalue corresponding to a+ is negative, perturbations leaving S = and 
changing only the magnitudes of the order parameters will decay exponentially in time. The 
eigenvalue corresponding to a_ is zero, so perturbations leaving S = and changing only 
the angles of the order parameters will be neutrally stable. 
Because £ = c + is, where 



' COs('0i) ' 

cos(V> 2 ) 
I cos(^m) J 



and s 



' sin(^i) ' 
sin(V> 2 ) 

^sin(^Af) J 



(24) 



the real vectors u satisfying = are exactly the vectors orthogonal to both c and s. For 
M > 3 this condition represents two constraints on u, so the space of vectors u is M — 2 
dimensional. Therefore there are M—2 independent eigenvectors of matrix A with eigenvalue 
0, and M — 2 with eigenvalue —(K(l + b) — 2). For M = 2, we have that tp2 = ipi + 7r. In 
this case the real and imaginary parts of t^u = give the same constraint, so the space of 
vectors u is one dimensional. 

Since A is a 2M x 2M matrix, our analysis accounts for all but four of the 2M eigenvalues 
of A. The previously identified eigenvectors span the space orthogonal to 





and 



s 



(25) 
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The matrix A must leave the space spanned by these four vectors invariant, because, for any 
M dimensional vector v, ££* T v and £*qv lie within that space. We can thus consider the 
restriction of A to that space, obtained by the transformation V~ 1 AV, where the 2M x 4 
matrix V is 

(c c s s \ 
- _ J ( 26 ) 
s —s c —c I 

and V~ x is a left inverse of V such that V~ l V is the 4x4 identity matrix. V~ 1 AV is a 4 x 4 
matrix whose eigenvalues are the eigenvalues of A that have corresponding eigenvectors 
in the restricted space. Using Mathematica, we find that these four eigenvalues are the 
solutions to the two quadratic equations, 



where 



A 2 + B ± X + C± = 0, (27) 



Kb 

B ± = (K(l + b) - 2) + ^°(M ± Pr 2 ), (28) 



Kb ( K \ 

C ± = — (1 + r 2 )(M ± P) ^-(1 + b) - lj 



(29) 



P= ^(^c-^^ + A^s) 2 



M M 



v EE cos [ 2 ^--^)]- ( 3 °) 

\ i=l i=i 



Note that < P < M. For all four roots of these quadratics to be real and negative, we 
need both B± and C± to be positive, and B\ > 4C±. Since B± is the sum of two positive 
numbers whenever K > K c , it must be positive. Additionally, C± is the product of five 
positive terms whenever K > K c , so it must be positive. In Appendix |A] we show that 
Bl > 4C± as well. 

We have thus established that the four remaining eigenvalues of matrix A are negative 
for all valid choices of the system parameters. Note that these eigenvalues correspond to 
perturbations 5a that are linear combinations of the vectors s and c. Such perturbations 
have 5S ^ 0. Since all four of these eigenvalues are negative, any perturbations that change 
the value of S will die out exponentially in time. 

This completes our analysis of the stability of the S = equilibria in which all order 
parameters have length r$. The equilibria are stable with respect to an M + 2 dimensional 
space of perturbations and neutrally stable with respect to an M — 2 dimensional space of 
perturbations. The neutrally stable perturbations all maintain the condition S = 0, so this 
space of equilibria is a stable set of solutions. 
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2. Equilibria with S = and one or more incoherent groups 



Suppose, in an equilibrium with S = 0, that m of the groups (denoted a — 1, 2, m) are 
incoherent, so that a\ = a<i = ■ ■ ■ = a m — 0, and the rest of the groups are coherent, with 
|a m +i| = |a m +2| = • • ■ = \olm\ — r o an d X^=m+i a °" = ^- We first consider the stability of 
this equilibrium for the case where m > 2. We perturb the equilibrium in such a way that 
oi\ = —Q2 (where a\ and a2 are small but non-zero) and all other a a are unchanged. This 



perturbation leaves S = 0. Linearizing Eq. (11), we find that 

da i 2 ( K 



it V2 a + 6)-lja.,, (31) 

For K > K c , this yields exponential growth for a\ and «2 [preserving the initial requirement 
that ai(t) = — a 2 (^)] and hence instability for the equilibrium. 

It remains to consider the stability of S = equilibria in which only one group is inco- 
herent. This is discussed in Appendix [Bj 

3. Equilibria with S ^ 

Since all the order parameters, as well as S, can be assumed to lie on the real line for 
equilibria with S ^ 0, we have S* = S and r* = r a for all a. We will consider strictly 
imaginary perturbations 5r a to these equilibria and show that they are unstable given those 
perturbations. Then if we define 6S = J2t=i ^ r c5 we a l so have 5S* = —5S and 5r* a = —5r a . 



Considering only the first-order differential terms of Eq. (11) and incorporating the above 
assumptions, we obtain 



d5r a 
dt 



l + ?-((l + b)(rl-l)-2r <7 bS) 



Kb 

Sr a -—(l + rl)5S. (32) 



Taking into account the equilibrium condition of Eq. (17), this equation becomes 



Kb dt y r. 
Now we rewrite this as a vector equation 

2 dfrr 



1 fM,VT 'r c + -) [S5r a - rjB\ . (33) 



(D + D-^iSt -lD)6r, (34) 



Kb dt 

where 1 is the M xM matrix with all entries equal to one, and D is the diagonal matrix D 



diagfrx, r 2 , r M }. We also define the matrices E = diagf^/ 1 + r 1 2 , a/1 + r 2 2 , J 1 + r 
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and W = S(D + D' 1 ) - E 2 D\D. Equation ([34} can then be written as 

2 d5a 



Kb dt 



WSa. 



(35) 



If we assume an exponential solution for Sr, then we see from Eq. (35) that the equilibrium 



will be unstable if the matrix W has any positive eigenvalues. The eigenvalues of W are the 
same as those of the matrix W = E~ 1 WE, which can be written as 



W = S(D + IT 1 ) - (ED)t(ED). 



(36) 



The matrix W is symmetric and has real entries, so the largest eigenvalue A max of W must 



satisfy 



A 



u T Wu 



mux _ ->- r 
Ill 



U x U 



(37) 



for any real vector u. We choose a vector u such that t(ED)u = 0. This condition can be 
written as 



M 



E"oV / ^ + ^ = 0, 



(38) 



<T=1 



where Uq is the <rth component of the vector u Q . This represents a one-dimensional constraint 
on u , so there exists an M — 1 dimensional subspace of valid Uq vectors. It now remains 



to show that we can find a uq for which the right hand side of Eq. (37) is positive. To 



accomplish this, we note that for any cr e {1, 2, M}, Eq. (17) can be solved for S, giving 

r CT (l + b) /rg - r 2 a 



S 



(39) 



where rg is defined as in Eq. (15). In order to satisfy the assumption that S > 0, we must 
have that r 2 . < rg if r a > and r 2 > rg if r CT < 0. This implies that the magnitude of the 



negative root of Eq. (17) is always greater than either of the two positive roots. Thus we 



cannot have S > unless at least two of the values of r a are positive; if all but one of the 
r CT values were negative, the sum of the r CT values would have to be negative. 

Let r ai and r U2 be any two positive order parameters. Choose the vector u so that all 



of the components Uq are zero when a ^ <r\ or 02, and choose Mq 1 and Uq 2 so that Eq. (38) 



is satisfied. Then Eq. (37) gives 



Ajrtax — 



u, 



<ri\2, 



ro 1 + r CTl x 



hi 



o 2 ) 2 ( 



r ao + r 



K x ) 2 + K 2 ) 2 
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(40) 



which is positive. Thus, this choice for uq shows that X m ax > 0. The matrix W therefore has 
a positive eigenvalue, implying that there exist perturbations 5a that grow exponentially 
in time. Consequently, any member of this class of equilibria (those with S 7^ 0) is unsta- 
ble. Thus the above analysis confirms our numerical observation that, among all possible 
equilibria, the only one that is not unstable is the one with all order parameter magnitudes 
equal to r and S = 0. 



IV. NONIDENTICAL GROUPS 
A. Formulation 

In order to make our model more realistic and reconcile the low-dimensional analysis with 



numerical simulations of the full system reported in Sec. HID we introduce small deviations 
in the frequency distribution centers A CT and widths u a , as well as in the coupling strengths 
K aa r of the oscillator groups. We let e be a small positive number and let u a = tuj a be 
the average natural frequency of the oscillators in the ath group. We also let e5 a be the 
deviation in the width of the frequency distribution of group a, and ek aa > be the deviation 
in coupling strength between groups o and a' . Thus 

A CT = 1 + e5 a (41) 

and 

{K + ekaa 1 : a = a' 
(42) 
-Kb + ek aa , :a^a' 

Letting r = et be a long characteristic timescale for the system, we then assume a configu- 
ration in which S = and a ff (t, r) = r e^ CT ^(l + ea£\t, r)), where ea^\t, r) represents a 
small, "fast" perturbation in a a , and the angles tp a , by virtue of their assumed dependence 



on r, are allowed to change slowly in time. Substituting this ansatz for a a into Eq. ( 11 ) and 
collecting terms that are first order in e, we get 

d(ai l) ) , m , K 



dt 



+ ai 1 ' + 



(1 + b) [r 2 aP* + (2r 2 - ljaW] - br e^S^> + -e~^S^ 

- i - «,,) - I f; Aw(r V^-^ - e -^~^)), (43) 

^ ' cr' = l 
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where = Y^=i r o e ^ c ' a< ^ ■ Note the similarity to Eq. (18). Just as in our stability 
analysis in Sec. |III E 1 , we can combine this system of coupled ODEs and their conjugates 
into a vector equation: 




(44) 



where A is defined as in Eq. (22), 



a 



(i) 



and the components of Q are given by 



,'dlpa 
-On- — I l — U r 



dr 



(1) 



M 



(45) 



(46) 



cr' = l 



We can find constraints on the value of Q by multiplying Eq. (44 ) on the left by ( u T 
where u is a vector satisfying = as in Sec. 



Ill El 



-u 



T 



The vector ( u T 



—if ) is a left 



eigenvector of A with eigenvalue 0, so the second term in Eq. (44) vanishes, leaving 



u 



T 




U 




(47) 



The right-hand side of Eq. (47) is effectively constant on the t timescale. If it is non-zero, 



then will diverge linearly with t on the fast timescale. We prevent this by demanding 
that 

' Q 

Q* 



u 



—U 



0. 



(48) 



or u^ilmQ) = 0, where ImQ denotes the imaginary part of Q. This equation holds for all 
valid choices of u if and only if Im Q is in the space spanned by the vectors c and s (see Eq. 



(24)). Thus, at each time r there are scalars X and Y such that 



dr 



1 M 

- ^ [fc<T<7'( r o + l)sin(^<7 - ^<j')] + X cosip a + Y smip a . 



(49) 



<t'=1 
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We find the coefficients X and Y by recalling that S = and therefore 

M M 

cos ipa = sin ^ CT = 0. 

a=l (7=1 

Taking the r derivative of this equation gives 

M A I M A I 

smip a -— = ^cost/v 



(7=1 



a=l 



dr 



0. 



(50) 



(51) 



Substituting Eq. (49) into Eq. (51 ), we obtain a system of two linear equations that uniquely 



specify X and Y at each time r: 



J2a=l COS Ipa) X + ( Yj a =l sin ^<J cos tye ) Y 



YLc=\ COS Ipu + \ Yla=l Yjtf=l k ^'{ r l + *) Sm(^ (J - tj) ff ,) COS V> CT 



1 v^M v^JU 



J2a=i sin cos ) X + ( J^Lj sin 2 if) a ) K 



- Eili ^ sin + \ Yfa=\ E"=i ( r o + !) singer - Vv) sin 



1 v^M 



(52) 



Solving Eq. (52 ) gives X and K as functions of the angles ip a . Insertion of these functions into 



Eq. (49) then gives the desired slow timescale evolution equation for if). Note that, because 



the values of 8 a are real, they do not appear in Eq. (52) and hence do not contribute to the 



0(e) slow timescale evolution. Also note that ^cos^V and J^sim/v are constants of the 



motion for this system (see Eq. (51)). 



The two linear algebraic equations in Eq. (52) are insoluble when the associated deter- 
minant, 



T = 



^cos 2 ?/^ ^^sin 2 Vv^ - sin ip a cos ip a 



{c 1 cjis 1 s) - (c 1 s 



(53) 



is equal to zero. This happens when c and s (the real and imaginary parts of £) are parallel, 
implying tan^i = tan ip2 = ■ ■ ■ = teaiipM- Thus T = occurs when all the if>„ are at two 
angles separated by tt (i.e., the order parameters lie on a single line through the origin). 
Because we still need r a ro and S 0, T = can only occur when M is an even number. 
For M = 2, T is always zero (if>i = ^2+vr). For M = 4 the four order parameters occur as two 



7r-separated pairs, so that T = whenever the two pairs coincide (Sec. IV B). For M even, 
as M increases beyond 4, it becomes unlikely that T = will ever occur. Our short/fast 
timescale expansion assumes that T = O(e ). Thus, in situations where T approaches zero, 



our expansion breaks down (when T becomes 0(e)). In Sec. IV C we will discuss the possible 
occurrence and implications of T becoming small in the context of numerical experiments. 
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B. The Examples of M = 3 and M = 4 



When M = 3, the space of vectors u satisfying = is one dimensional, spanned by 
the vector 



M 



M 

1 

W 



(54) 



In this case, Eq. (48) gives 



£ 

<r=l 



1 3 
+ ^(^o + !) ^ CTCT ' sin (^ ~ ^ CT ')] 



o-' = l 



0. 



(55) 



Because we must have 5 = at each time step, the phases if) ff must be separated by 27r/3, 
and all three phases must change at the same rate. If we assume ip2 — ^i — ^3 — ^2 — 27r/3, 



Eq. (55) becomes 



<9Vv 



W w + tUfc, 



where 



v3, 2 

\r 

3 12 



Wi + u 2 + 003 



r\ + 1) [k 12 + A; 2 3 + &3i - &21 - k. 



32 — tl3l 



(56) 



(57) 



and a = 1,2, or 3. We thus see that the phases ip a change at the same constant rate given 



by the right-hand side of Eq. (56). Note that, if the perturbed coupling matrix is symmetric, 



k a(T i = k a i a , then u> fc = 0, and the rotation is at the average frequency w w . 

In the M = 4 case, S = is obtained by having two pairs of order parameters whose 
phases are separated by ir. We therefore assume without loss of generality that ^3 = + rr 
and ipA = ip2 + We can then express the vectors u independently of the angles ip a . One 
possible choice is 



Ui 




1 

w 



and m 2 



1 




(58) 



It is clear that these vectors satisfy £ u — given the above conditions on the phases ip a . 



If we use vector u\ in Eq. (48), we obtain 

dip 



^- Wl + I^ lCT (r 2 + l)sin(^^ (59) 



(7=1 



(7=1 
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In order to preserve S = 0, it is necessary that ipi an d ip3 change at the same rate. Letting 



dip 3 /dr = dipi/dr in Eq. (59) and simplifying yields 

dlpi 



dr 



+ -(rg + l)(fci2 - fc 32 - fc M + ^34) sin(V'i - 2 ) 



(60) 



Using vector «2 in Eq. (48) gives a similar ODE for ip 2 : 
dip 2 



dr 



U2 + 2 UA + -^{tq + l)(k 2 i - k 4 i - k 23 + k 43 ) sin(^ 2 



0i 



0. 



(61) 



If we subtract Eq. (61) from Eq. (60) and define Aip = ip\ — ip 2 , we obtain an ODE for the 
evolution of Aip: 



dAip 



-tt + UsmAip = 0, 



(62) 



where fl = {uj\ — W2+W3 — u±)/2 and H 



l){k 12 + k 2 i + k u + h 3 - ku - kn -k 23 -k 32 )/4:. 



Once Aip = (vpi —ip 2 ) is found from Eq. (62), ip\ and ijj 2 can be determined by inserting Aijj 



into (60) and (61). Note that f2 depends only on the frequency perturbations co a , while the 



term H responsible for interaction between groups depends only on the off-diagonal coupling 



perturbations, k aa i for o ^ a'. The solution to Eq. (62) depends on the relative sizes of 



and H. If > the solution is periodic, and the time derivative of Aip always has the 
same sign. If < \7i\, then Aip eventually reaches a constant value, at which point all 
groups rotate at the same speed. 

In this analysis of the M = 4 case we assume that the pairs remain the same, i.e., that 
groups 1 and 3 always remain paired. However, it is possible for the groups to switch 
partners when the group order parameters overlap (A-0 = or 71). In this case the vectors 
u\ and u 2 will change to reflect the new pairing. This phenomenon will be discussed in Sec. 

wu 

C. Numerical Results 



Returning to our numerical simulation of the low-dimensional system, we introduced the 
small perturbations co a , 5 a , and k acT i, and integrated Eq. (m). When these perturbations were 
included, the system no longer approached a steady-state equilibrium. Rather, the group 



order parameters behaved as in the high-dimensional simulations of Sec. |III D[ establishing 
a configuration in which S ~ and then evolving in the complex a-plane in a way that 
kept S near zero and \a a \ nearly r at all times. The similarity to the high-dimensional 
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system indicates that the behavior of the full system can be well understood by studying 
the low-dimensional case with nonidentical groups. 



Figure [4] shows a numerical example for a case with three groups (M = 3). In this 
example the parameter perturbations from exact symmetry are 



cojx = 0.0158, eu 2 = 0.0060, eco 3 = -0.0033, 



e<5i = .0252, e5 2 = -0.04 32, e5 3 = -0.018* 



(63) 



ek 



0.0321 0.0034 0.0148 
0.0270 -0.0224 0.0180 
-0.0090 0.0080 -0.0306 



The initial conditions (Fig. |4]^a)) are such that S is not initially zero. However, as time 
increases, S is observed to move quickly towards zero (Fig. g|b)), but due to the O(e) 
parameter perturbations from symmetry, its magnitude, although small, does not become 



exactly zero. Also, the order parameter magnitudes r a approach r given by Eq. (15) (Fig 



||c)), but with some O(e) deviation. Figure |4](d) shows the time evolution of the order pa- 



rameter phases ipi,ip2, and ip3, which, as predicted in Sec. IV B approach an approximately 
evenly spaced, uniformly rotating configuration. Figure |4](e) shows the time evolution of the 
average rotation rate dip/dr, where ip = + ip2 + V ; 3)/3, showing that dip/dr approaches 



the theoretical prediction + Wk given in Eq. (56). 
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FIG. 4. Results for a case with three groups {M — 3), obtained by integrating Eq. ^ with w ff ,A„, and 



ifo-CT' defined as in Sec. IV A and the perturbations aj ai S a , and k aa i given by (63). (a) The initial conditions, 
ai(0) (red dot), 0:2(0) (green dot), 03(0) (blue dot); (b) the evolution of \S\\ (c) evolution of r CT with respect 
to tq; (d) evolution of ipi, and ^3 (in radians); (e) evolution of di/i/dr, where = (?/>i + ip2 + V'3)/3. In 



(e), the horizontal line denotes the value w — u> w + Wk given in Eq. (57 1 



When M = 4, numerical simulations exhibited both types of behavior predicted by Eq. 



(62): when the values of cu a and k aa i were such that > \H\, Atp evolved periodically, 
and when |0| < [H\, Aip approached a fixed value. However, the behavior was often more 
complicated than this, due to a pair switching phenomenon. If, for example, the initial 
configuration was such that 03 = ipi + it and -04 = 02 + 7T, it was sometimes possible for 
the groups to switch partners midway through the simulation (see Fig. ga) at t « 130 and 
t w 860), so that ipi became paired with -02 (^2 — V'l + 7r ) an d ^3 became paired with -04 
(■03 = ^4 + 7r). In Fig. ga), for example, when t < 130, the two pairings are red with black 
and blue with green; when 130 < t < 860, the pairings are red with green and black with 
blue; and when t > 860 the pairings are red with blue and black with green. The values of 
Q and % depend on the pairings of the groups, so it is possible to have \Q\ > \H\ before a 
pair switch and \ < \H\ afterwards. This is what happens at t ~ 860 in Fig. ga). Before 
t rj 860, |0| > \H\. After t ft* 860, [H\ > and, for t > 860, Aip consequently approaches 
a fixed value, after which the whole configuration rotates rigidly. Both the value of A-0 and 
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the value of the rotation rate after t ~ 860 are correctly predicted by the fixed point of Eq. 



1500 
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(b) 
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1000 



(e) 



(f) 



500 



FIG. 5. Numerical results for three cases with four groups (M = 4). (a, c, e) The evolution of ipi, ip2, ip3, 
and "04 (in radians); (b, d, f) the corresponding evolution of T. The pair switches in each simulation are 
possible only when the value of T nears zero. In (a) and (b), this occurs at t « 130 and t « 860. In (c) 
and (d), the pairs switch in a periodic manner throughout the simulation. In (e) and (f), we see that T ~ 
several times in the simulation but only one pair switch occurs, at i ~ 35. 



The pair switching phenomenon is related to the failure of Eq. (52) to completely de 



termine the system's behavior. Figure [51(b) shows the value of the quantity T defined in 



Eq. (53) as a function of time. We observe that T each time a pair switch occurs. 



(For t > 860 in Fig. [5|a, b) the rigid rotation of the configurations implies that T becomes 
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constant and that no pair switchings can occur after t ~ 860.) 

Another M = 4 example is shown in Figs. [5]^c, d). In this case we observe that the 
switching assumes an apparently persistent time-periodic pattern. For example, considering 
the group plotted in green in Fig. |5|c), we see that it is initially paired with red, then 
switches to being paired with blue, then with black, then back to red, etc. Correspondingly, 
Fig. |5](d) shows that T becomes periodic in time. We observe that the periodic patterns in 
Figs. |5](c, d) continue for the length of very long numerical runs. 

Note that it is also possible to have T w when the two pairs pass through one another 
without switching. For example, in Fig. [5|e, f), T approaches zero many times throughout 
the simulation, but only the first of those occurrences (t w 35) corresponds to a pair switch. 

Figure [6] shows results for a numerical simulation with M — 5. When M — 5, it is 
impossible to have T = as long as the system maintains a configuration in which S ~ 
and r a m r for all a. Thus, the evolution of the phases ipo- (Fig- [6^ a )) always obeys Eq. 
(49). Note that the value of T (Fig. [6tb)) is positive throughout the simulation. 

The results of a simulation with M = 6 can be seen in Fig. [7| When M = 6 it is possible, 
though unlikely, for T = to occur. We observe that T decreases significantly many times 
throughout the simulation, but it never approaches zero as closely as in the M = 4 case. This 



means that Eq. (52) remains solvable for the duration of the simulation, so the dynamics of 



the phases can be described by Eq. (49). However, T may sporadically become quite small 



in some M = 6 simulations, so Eq. (49) is not guaranteed to hold in all M — 6 cases. 



Because of the difficulty in expressing equations of motion akin to Eqs. (56) and (62) 
for M = 5 and 6, it is not straightforward to characterize the possible behaviors of these 
systems. For example, in some M = 5 and M = 6 simulations, the system exhibited periodic 
behavior. In other simulations, the order parameters became locked in a static configuration, 
similar to the one seen in Fig. |Ja) after t w 860. 



V. CONCLUSION 



We have considered the dynamics of M interacting groups of coupled oscillators, with 
"attractive" coupling within groups and "repulsive" coupling between groups. If the number 
of oscillators in each group is large, this system can be described using the low dimensional 
ansatz of Ott and AntonseiP, which reduces the problem to a set of M ODEs for the M 
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1000 2000 3000 4000 5000 

Time 

FIG. 6. Numerical results for a case with five groups (M — 5). (a) The evolution of ipi, -02, "03, "04: an d "05 
(in radians); (b) the corresponding evolution of T. The value of T remains large throughout the simulation, 
in contrast to the cases with M = 4 in Fig. [5| 

complex group order parameters. Based on our analyses the typical dynamical evolution for 
such a system is as follows. There is an initial phase of evolution in which members of the 
same group, which are attracted to each other, become synchronized and the group acquires 
a nonzero, complex order parameter. In the symmetric case, where the group properties and 
interactions are all the same, the group order parameters relax on the same time scale to a 
manifold of neutrally stable equilibria in which all of the group order parameters have the 
same magnitude, and due to the repulsion between groups, the sum of the complex order 
parameters is zero. Other equilibria also exist in this symmetric case. However, we have 
shown them to be unstable. 

In order to gain insight into the general non-symmetric case, we then introduced small 
perturbations to the group frequency distributions and coupling strengths. The effect of 
these small asymmetries is to introduce a slow time evolution to the phases of the group 
order parameters of the neutrally stable equilibria found in the symmetric case. To study 
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Time 



FIG. 7. Numerical results for a case with six groups (M = 6). (a) The evolution of ipi, ip2, ip3, ?p4, ip5 and 
ipQ (in radians); (b) the corresponding evolution of T. 

this we used a slow/fast time scale analysis and obtained ODEs describing the slow-timescale 
motion of the order parameter phases. In the M = 3 case we observed that the group order 
parameters remained evenly spaced and rotated at a constant rate in the complex plane. In 
the M = 4 case the order parameters formed two pairs where members of the pairs differed 
in phase by tt. Depending on the asymmetry parameters the relative phases of the pairs 
either became locked or evolved periodically in time. In the M = 5 and higher cases the 
behavior of the order parameter phases became more complicated, showing a mixture of 
irregular and periodic motion. Our conclusions above are supported by both analysis and 
numerical simulations of the low dimensional and high dimensional equations. 
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Appendix A: Demonstration that B± > 4C± 

To determine whether B± > AC±, we take its derivative with respect to P to obtain 

d{Bl ap AC±) = ±KV ° [-(l + b) + l(M± Pr 2 )] . (Al) 
Note that the second derivative with respect to P is always positive, so all extrema of this 



function are minima. If we choose the "— " option in Eq. (Al), we find that this derivative 
is zero when 

p = Tq 2 (M - 2(1 + b) fb). (A2) 
Then there are two cases depending on the value of M. If M > 2(1 + b)/b, then the value 



of B 2 _ — 4C„ is positive at the minimum, as hoped. If M < 2(1 + b)/b, then Eq. (A2) gives 
a negative value, but P must be at least zero so the minimum value of B 2 _ — 4C_ occurs 
outside the possible range of values for P. Therefore, B 2 _ — 4C_ is increasing on the interval 
P G [0, M] so we need only check its value at P = to ensure it is always positive. Its value 
at P = can be expressed as ((K(l + b) - 2) - KMb/2) 2 + KMb{K{\ + b) - 2)(1 - r 2 ), 
which is positive. 



Similarly, if we choose the "+" option in Eq. (Al), then the derivative is zero when 

P = r 2 {2(1 + b)/b- M). (A3) 
Again there are two cases depending on the value of M. If M > (1 + b)/b, then the 



value of B 2 + — 4C+ is nonnegative. If M < (1 + b)/b, then Eq. (A3) gives a value greater 
than M, so the minimum value for B\ — AC + again occurs outside the possible range of 
values for P. Therefore, B 2 ^ — 4C + is decreasing on the interval P G [0, M] so we need 
only check its value at P = M to ensure it is always positive. Its value at P = M is 
((A(l + b) - 2) - KMb(\ + r 2 )/2) 2 which is clearly positive as well. Thus B\ > 4C±. 

Appendix B: Instability of S = equilibria with exactly one incoherent group 

Numerical simulations imply that equilibria in which S = and exactly one group is 
incoherent are unstable. However, because the stability analysis of this case is algebraically 
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intensive and the numerical evidence for instability is compelling, we prove instability only 
for the more tractable special case in which the nonzero order parameters are evenly spaced 
on a circle of radius r in the complex plane. Defining A8 = 2ir/ (M — 1) and letting 5a a be 
small complex numbers, we suppose that 

a a = e i(rA0 (r o + 8a a ) (Bl) 



for a = 1,2,...,M — 1, and au = 8(Xm- Substituting these quantities into Eq. (11) and 
collecting terms that are first order in 5a a gives 

^ = -f rl{l + b){8a a + Soft ~ f & {8Se*>*° - rg^e***) (B2) 



for a = 1,2,...,M-1, and 



dt 2 0V ; 2 



(B3) 



In order to describe the collective behavior of the perturbed system, we define 

M-l M-l 



e iaA0 8a a and a* = e iaAe 8a* a , (B4) 



a = 

(7=1 (7=1 



and multiply Eq. (B2) on both sides by exp(?crA#) and sum over all values of a from 1 to 



M — 1. Assuming d(a)/dt = Act for some scalar A, we have 



and 



where 



2 

We also assume d(8aM)/dt = \5chm, so that 



(A + A )ct = -A a* - —bM5S (B5) 



(A + A )ct* = -A ct - —br 2 M5S, (B6) 



A = ^r 2 (l + 6). (B7) 



2 

Because is the sum of the perturbations, we have 



(A - X )6a M = -^-SS. (B8) 



5S = a + 5o>m- (B9) 
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Substituting this into Eq. (B8) and rearranging terms, we get 

A -Ac 



a + 5a 



M 



a 



A - A + Kb/2 



(BIO) 



Replacing 5S with (BIO) in Eqs. (B5) and (B6) yields 



and 



(A + A ) + 



(A + Ao)a* 



KbM 



A-Ar 



2 A - A + Kb/2 



a 



-A a* 



-A + 



KbM 



A-A 



2 A - A + Kb/2 
Combining these two equations and collecting powers of A gives 



a. 



(Bll) 



(B12) 



Kb 



A +A —(M+l) + X )+X[Kb\o + —Ma 2 Xo 



Kb 



2X1 



Kb 



MX 2 (l + a 2 ) = (B13) 



The value of this polynomial at A = is — {KbM /2)X 2 ) (l-\- ofy . This is a negative number 
when K > K c , because and Ao are positive. Because the polynomial in Eq. (B13) is 
negative at A = and its value tends to positive infinity as A increases, the polynomial must 
have a root with A > 0. The presence of a positive root indicates exponential growth of the 
perturbations, which means that this equilibrium is unstable. 
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